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GENERAL RELATIVITY AND SATELLITE ORBITS 

by 

David Parry Rubincam 
ABSTRACT 

The general relativistic correction to the position of a satellite is found by 
retaining Newtonian physics for an observer on the satellite and introducing a 
r~ 3 potential. The potential is expanded in terms of the Keplerian elements of 
the orbit and substituted in Lagrange r s equations. Integration of the equations 
shows that a typical earth satellite with small orbital eccentricity is displaced 
by about 17 cm from its unperturbed position after a single orbit, while the 
periodic displacement over the orbit reaches a maximum of about 3 cm. The 
moon is displaced by about the same amounts. Application of the equations to 
Mercury gives a total displacement of about 58 km after one orbit and a maximum 
periodic displacement of about 12 km. 


iii 



CONTENTS 

Page 

INTRODUCTION 1 

DERIVATION OF THE EQUATIONS OF MOTION 2 

NEWTONIAN FORMULATION 3 

EQUATIONS OF CELESTIAL MECHANICS 4 

INTEGRATION OF THE EQUATIONS 9 

THE ORBIT AS SEEN FROM THE GROUND 14 

ALTERNATIVE VIEWPOINT 19' 

CONCLUSION 20 

ACKNOWLEDGMENTS 20 

REFERENCES 21 

APPENDIX 22 

LIST OF ILLUSTRATIONS 

Figure 1. Orientation of the Orbit 11 

Figure 2. The total displacement of Beacon Explorer C due to the 

relativistic potential. The unperturbed position is at point 
O. The perturbed position is at point A after one revolution 
of the unperturbed orbit. The numbers in the diagram refer 
to centimeters. The diagram itself is one -half actual 

size 12 

Figure 3. Total displacement of the moon over one orbit. The 

numbers in the diagrams refer to centimeters. One-half 

actual size 13 

Figure 4. Total displacement of Mercury over one orbit. The distorted 
shape of the spiral is due to the large eccentricity of the 
- orbit. The numbers in the diagram refer to kilometers ... 15 

Figure 5. Secular displacement of Mercury. The numbers in the 

diagram refer to kilometers 16 



ILLUSTRATIONS (Cont.) 


Page 

Figure 6. The periodic displacement of Beacon Explorer C (top) 

and the moon (bottom). The sense of rotation is counter- 
clockwise. The numbers refer to centimeters. Actual 

size. 17 

Figure 7. The periodic displacement of Mercury keeping powers of 
e < 2 (top) and £ 1 (bottom). The top diagram gives a 
more accurate picture of the periodic displacement than the . 
bottom diagram. The numbers refer to kilometers 18 

LIST OF TABLES 

Table 1. Eccentricity functions. From Kaula (1966), Caputo (1967), 

and Cayley (1861) 5 

Table 2. Secular rates of change of the argument of perigee co and 
mean anomaly M for Beacon Explorer C, the moon, and 

Mercury 8 

Table 3. Summary of displacement data for Beacon Explorer C, the 

moon, and Mercury 9 


vi 



GENERAL RELATIVITY AND SATELLITE ORBITS 


INTRODUCTION 

The primary purpose of this work is to investigate the effect of general relativity 
on the orbits of artificial satellites; but the results may be applied to any body 
of negligible mass orbiting about a massive, spherically symmetric object. In 
particular, we will discuss the moon orbiting around the earth and the planet 
Mercury orbiting around the sun. 

Past attacks on the problem have centered around solving the equation (see 
Ghaffari, 1970 and references contained therein): 


d 2 u GM 3GM 2 

+ u = + u z 

dc£ 2 h 2 c 2 

Our technique will be to back off from this equation a little to a point where we 
may interpret the equations of motion in the following way: the geometry of 
space is Euclidian and the physics is Newtonian. The price we pay for this 
approach is that we must modify the law of gravity and introduce an extra 
(relativistic) potential. This poses no particular problems, however, since the 
relativistic potential now becomes a disturbing function susceptible to the 
methods of celestial mechanics. In particular, the potential may be expressed 
in terms of the Keplerian elements of the orbit and substituted in Lagrange r s 
equations. The equations may then be integrated to give the osculating elements 
of the orbit. 

The technique has a sound philosophical basis. Even though Einstein (and others) 
developed general relativity within the framework of non-Euclidian geometry, 
we can, however, obtain an equivalent description of the world by retaining 
Euclidian geometry and modifying the laws of phj^sics. This was discussed by 
Poincare'' (1905) and clearly explained by Carnap (1966). Convenience dictates 
the point of view we choose. Poincare felt mankind was so accustomed to 
Euclidian geometry that it might never abandon it in favor of non-Euclidian 
geometry, even though the latter point of view might represent a simpler picture 
of the world. Einstein and physicists in general, however, adopted the non- 
Euclidian approach for reasons of conceptual clarity and mathematical elegance. 
Indeed, it is doubtful general relativity could have been developed without it. But 
we wall follow Poincare' and introduce an extra Newtonian force, since it results 
in an elegant description of the motion of a satellite. 
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(Elegant for satellites but perhaps for not much else. For instance, if we measure 
the circumference and radius of a circle about the earth we discover that their 
ratio is not tt , Hence we would need some laws about the expansion and contraction 
of meter sticks. This particular problem is ignored here since the displacements 
we are concerned with are so small that this effect may be neglected.) 


DERIVATION OF THE EQUATIONS OF MOTION* 

Let us consider the motion of a body with negligible mass about a massive central 
object. We will call the two bodies satellite and earth, respectively, since we 
are primarily concerned with the motion of artificial satellites about the earth. 

The geometry of spacetime in the neighborhood of a spherically symmetric 
earth is given by the Schwarzschild line element (Tolman, eq. 82.9): 

ds 2 = - — - r 2 dd 2 - r 2 sin 2 6dcp 2 + (l - ?]c 2 dt 2 (1) 

V () 

Here ds is the interval of proper distance, (r, 6 , <p) are polar coordinates, t is 
the coordinate time, and M E is the mass of the earth. The speed of light c and 
universal constant of gravitation G are explicitly retained, in contrast to the usual 
procedure of setting G = 1 and c = 1. 

The satellite will follow a geodesic in the Schwarzschild geometry according to 
the geodesic equation (Tolman, eq. 83.1): 


d 2 x^ 

ds 2 


+ { )jlv , a} 


dx^ 

ds 


dx v 

ds 


= 0 


(2) 


where r - x 1 , £ - x 2 , 0 = x 3 , t = x 4 , and {jxv, a} is the Christoffel symbol. 

One may derive from equations (1) and (2) the equations of motion along with 
two constants of motion, k and h (Tolman, eqs. 83.10-11): 




^ 2 ] - 
c 2 dr 2 / 


(k 2 - 


1) c 2 


(3) 


★ 

Our treatment summarizes that of Tolman (1934). 
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(4) 


Here ds — cdr, so that dr is an element of proper time as measured by a clock on 
the satellite (Tolman, pg. 207). Angle 9 does not appear in equations (3) and (4) 
since Q has been set equal to V 2 without loss of generality. Hence the satellite 
remains in a fixed plane passing through the center of the earth. 


NEWTONIAN FORMULATION 

Substitution of (4) into (3) and dividing by 2 yields 


T + V N + = constant (5) 

where 



and 


gm e 

r 


V cr 


GM„h 2 
c 2 r 3 


(6) 

(7) 


Consider an observer on the satellite. His space coordinates are (r, 6 , <£>) and 
he measures time r with his clock. If the observer assumes his space is 
Euclidian and his physics is Newtonian, then T is the kinetic energy per unit 
mass oi the satellite and V N is the ordinary Newtonian potential. V GR is the 
general relativistic potential which we are now forced to introduce. 

Equation (5) now represents conservation of energy, while equation (4) represents 
conservation of angular momentum. (That angular momentum is conserved is 
easily seen from equations (6) and (7); both potentials represent central forces.) 
Hence, from the point of view of the satellite , it moves through Euclidian space 
under the action of the total potential V N + , with conserved energy and 

angular momentum. 
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equations of celestial mechanics 

With our Newtonian approach in hand, we are now ready to apply the methods of 
celestial mechanics. 

A satellite moving under the influence of V N only will describe an ellipse with 
constant orbital elements (except for M 0 ) a 0 , e 0 , i 0 , M 0 , , fi 0 . A satellite 

moving under the added influence of a disturbing function R will have osculating 
elements a, e, i, M* co, 0 which change in time according to Lagrange's equations 
(Brouwer and Clemence, 196 1; Blanco and McCuskey, 1961). 


where 


da 
dt ~ 

2 BR_ 
na BM 

de _ 
dt 

tL - e2 j/r 

na 2 e ^ 

di 

1 

dt 

na 2 /l - e 2 : 

dn _ 

1 

dt "■ 

na 2 /I - e 2 

do; 

dr “ 

✓ 1 - e 5 BR 
na 2 e 3e 

dM _ 
dt 

2 c)R 

na Ba 


BR 

BM 


BrI 

~ dojj 


1 fBR 

< COS 1 

■ 2 • - IBO 

- e s i n i ^ 


cos 1 


na' 


rva 2 e de 


br\ 

dco J 


BR 
B i 


/gm e 

n 

a 3/ 2 


For the case under discussion we set r = t and - R. 
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Our task now is to express V GR in terms of the orbital elements and substitute 
in Lagrange's equations. This may be elegantly done by noting that (Kaula, 1966 
Caputo, 1967): 


1 


l 


+ 1 


cos{(^ - 2p) (co 4 f) 4 m(fi - 8)} 



00 



Go (e)cos{( / t-2p)(4>4( / C'-2p4q)M 

'Vpq 

4 m(fl - (9)} 


Here f is the true anomaly, 8 is the Greenwich sidereal time, and the (e) 
are the eccentricity functions. Tables of G^ (e) may be found in Kaula (1966), 
Caputo (1967), and Cayley (1861); see also TaSle 1. 


Table 1 

Eccentricity functions. From Kaula (1966), Caputo (1967), 
and Cayley (1861). 



Setting = 2; p = 1 ; and m = 0, we obtain 


1 


1 


r 3 a 3 


/*~| G 21q (e) cos(qM) 

q = - co 


so that 


R = 


GM E h 2 


/ G 21q (e) cos(qM). 

q sr - 0D 
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The areal-velocity constant h may be evaluated in terms of the orbital elements. 
From considerations of the osculating ellipse we find (Blanco and McCuskey, 
pg. 133): 

h 2 ^ GM E a(l - e 2 ). 

Substitution of the disturbing function into Lagrange's equations yield 


di 

dt 


= 0 


( 8 ) 


dn 

d t 


rii- 0 


(9) 


(GM e ) 1/2 h 2 (1 _ e 2)l/2 
d t „2 


ea 


7/2 


GO 


(e) cos (qM) 


q = ~co 


(10) 


dM 

dt 


(GM r ) 1/2 h 2 


. 1/2 


1 T G' (e)cos(qM) 


q = _ oo 


) C 21 q(e) cos(qM) 


(GM e ) 1/2 


. 3/2 


( 11 ) 


2h 2 (GM ) 1/2 1 

= — > qG aiq (e) sin(qM) 

dt _2 ,-. 5/2 L 1 21q 


q = — 00 


( 12 ) 


de (GM E ) 1/2 h 2 t _ e 2 


dt 


ea 


7/2 


<l G 2i q ( e ) s in (<lM) 


~ _ 00 


( 13 ) 


The prime on G 21q (e) denotes differentiation with respect to e. 



We see from equations (8) and (9) that the inclination i and node fi remain constant. 

We obtain the secular rate of change of the elements by examining the terms for 
which q = 0: 



d^j 

_dt_ 


s 


(G!VL.) 1/2 h 2 (l - e 2 ) 1/2 


c 2 ea 7 ' 2 


G 2 io ( e ) 


(14) 



(GM e ) 1/2 h 2 

c 2 a 7/2 



^ 21 0 ( e > 


(1 - e 2 ) 




(15) 

S,N 


Here the subscripts S, GR, and N mean "secular", "general relativity" and 
"Newtonian", respectively. 

There is no secular change in the semimajor axis a or eccentricity e. 

The well-known expression for the rate of rotation of the argument of perigee qj 
may be obtained by noting that 


G 2io( e > = 0 - -O 


2 \” 3/2 


G W e > = 


3e 


(1 - e 2 ) s/2 

and h 2 = GM a(l - e 2 ). Substituting these expressions into equation (14), we get 


d co 


dt 


3(GM ) 3/2 


, 5/2 


(1 - e 2 ) 


This agrees with the usual expression obtained by other means (Tolman, 1934; 
Bergmann, 1942). Table 2 contains the secular rates of change of the argument 
of perigee for the satellite Beacon Explorer C, the moon, and the planet 
Mercury. 
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Table 2 

Secular rates of change of the argument of perigee co and mean anomaly M for 
Beacon Explorer C, the moon, and Mercury 


Orbiting body 

Central 

object 

gm e 

C 2 

a 

e 

day 

dt 

s 

[dM] 

Ws.GE 

Beacon Explorer C 

Earth 

0.443 cm 

7.502 x 10 8 cm 

0.025 

0.03 07’ '/day 

0.0307 "/day 

Moon 

Earth 

0.443 cm 

3.844 x lO 10 cm 

0.055 

1.63 x 10’ 6 ’’/day 

1.63 x 10 -6 ’’/day 

Mercury 

Sun 

1.477 km 

5.791 x 10 7 km 

0.206 

43 ’’/century 

42 '’/century 




From equation (15) we may write 

rdMi ( gm e > 1/2 h2 


~2„7/2 
S , GR c a 


6G 210 


(e) - G' 2I0 (e) (Llf! 


3(GM_) 1/2 h 2 


3(GM e ) 3/2 


c 2g7/2(l __ e 2^3/2 c 2 a^ /2 (l - e 2 ) 1/2 

Values for the secular rate of change of the mean anomaly M for the above men- 
tioned bodies may be found in Table 2. 

Table 3 

Summary of displacement data for Beacon Explorer C, the moon, and Mercury 


Orbiting body 

Displacement after 

Maximum periodic 

one orbit 

displacement 

Beacon Explorer C 

16.7 cm 

3.1 cm 

Moon 

16.8 cm 

3.1 cm 

Mercury 

58.1 km 

12.4 km 


INTEGRATION OF THE EQUATIONS 

If we substitute the elements of the unperturbed orbit a Q , e , i Q> M Q , into 
the right side of equations (10)-(13) and set 


(GM„) 1/2 


then the equations may be integrated to give 


/CM\ 3 /GM E \2(l-e 2 y 

^ + M o + 

\c 2 j (1 - e 2 ) a 0 lc 2 / e o a o 


sin qM 0 


V"— i s in i 

^G' 21 (e 0 )__ 
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as the approximate expressions for the elements of the perturbed orbit. The 
constants of integration have been adjusted so that the perturbed and unperturbed 
elements are identical at perigee. 

The displacement in the position of the satellite due to the relativistic potential 
may now be found. This was accomplished with the computer program given in 
the Appendix. The program takes the Keplerian elements of the perturbed orbit 
given by equations (16)-(19) and converts them to Cartesian coordinates and 
velocities x, y, z, x, y, z. The vector r = (x, y, z) gives the perturbed position 
of the satellite. The process is repeated for the unperturbed orbit, obtaining 

tVn a nnaitinn Tror-trir r — iv -xr >7. V Thf* Hi ff A 7’ — T* — V „ ETiveS the 
i *□ v^o’^O* “0 - - u 

displacement of the satellite due to the relativistic potential. 


The orbits are taken to lie in the xy plane as shown in Figure 1, so that z = z 0 = 0 
and z = z 0 =0. The x axis lies along the line of perigee of the unperturbed 
orbit, and the mean anomaly increases in the positive (counterclockwise) sense. 
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Figure 1. Orientation of the Orbit. 


The program uses double precision variables and can consider powers of e < 2 
or < 1 in the eccentricity functions and their derivatives, depending upon the 
choice of the programmer. 

Figure 2 plots A"r for a typical earth satellite (Beacon Explorer C). The un- 
perturbed position is at the origin (point O in the figure). The perturbed and 
unperturbed positions are coincident at perigee (point O) and the perturbed 
position moves away from the unperturbed position in a spiral as time progresses. 
After one revolution the perturbed position is at point A, about 16.7 cm from 
the unperturbed position. The displacement of the moon is given in Figure 3. 

Here the displacement is about 16.8 cm at the end of one revolution. 

Figure 4 gives A 7 for Mercury. The rather distorted shape of the spiral is due 
to the large eccentricity of the orbit. Mercury is displaced by about 58.1 km at 
the end of one revolution. 

The secular displacement of Mercury is shown in Figure 5. Here the periodic 
terms have been omitted. 
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Figure 2. The tota] displacement of Beacon Explorer C due to the relativistic potential. The unperturbed position is at 
point O. The perturbed position is at point A after one revolution of the unperturbed orbit. The numbers in the diagram 
refer to centimeters. The diagram itself is one-half actual size. 
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Figure 3. Total displacement of the moon over one orbit. The numbers in the diagram refer to centimeters. One-half actual size. 
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The periodic displacement for Beacon Exploer C (top) and the moon (bottom) 
are given in Figure 6. In each case the maximum displacement (maximum |Ar |) 
is about 3.1 cm. The sense of rotation is counterclockwise. 

The periodic displacement for Mercury is shown in Figure 7 (top). The maximum 
distance from the unperturbed position is about 12.4 km. 

A summary of the numerical data for Figure 2-7 is given in Table 3. 

Powers of e < 2 were retained in the eccentricity functions and their derivatives 
for the computations for Figures 2, 3, 4, 6, and 7 (top). The periodic displacement 
for Mercury with powers of e < 1 retained is given in Figure 7 (bottom). The 
great difference in the shapes of the curves emphasizes the importance of keep- 
ing many terms in the eccentricity functions and their derivatives when the 
eccentricity is large. 

It is interesting to note that the displacements are the same for any two bodies 
orbiting the same massive central object, regardless of the values of the semi- 
major axes, so long as the eccentricities are the same. This is apparently due 
to the appearance of a 0 in the denominators of the corrections to co 0> M 0 , and e 0 
(equations 16, 17, and 19), As a 0 changes scale, co Q , M 0 , and e 0 change in such 
a manner as to compensate for it. This explain why the displacement of 
Beacon Explorer C and the moon are very nearly the same. 


THE ORBIT AS SEEN FROM THE GROUND 


So far we have proceeded from the point of view of an observer on the satellite. 
An observer on the ground sees a somewhat more complicated set of forces 
acting than does the observer on the satellite; this may be verified by examining 
the Schwarzschild metric from the point of view of the ground observer. We will 
not pursue this very far, except to say that both observers will agree on the 
track of the satellite across the sky. That this is so may be seen by imaging 
space to be laced with coordinate lines. The two observers must agree that the 
satellite arrives at various points in the coordinate system as the satellite moves 
around the earth. This line of reasoning is implicit in the equation given in the 
Introduction; the solution of the equation gives r (= u" 1 ) as a function of 0, which 
is the same regardless of who is looking at the satellite. Hence when both the 
observer on the ground and the satellite plot out the orbit, they will get the 
same answer. 
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Figure 4. Total displacement of Mercury over one orbit. The distorted shape of the spiral is due to the large eccentricity of 

the orbit. The numbers in the diagram refer to kilometers. 
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Figure 5. Secular displacement of Mercury. The numbers in the diagram refer to kilometers. 
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Figure 6. The periodic displacement of Beacon Explorer C (top) and the moon (bottom). The sense of rotation is 
counterclockwise. The numbers refer to centimeters. Actual size. 
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Figure 7. The periodic displacement of Mercury keeping powers of e ^ 2 (top) and (bottom). The top diagram gives a 
more accurate picture of the periodic displacement than the bottom diagram. The numbers refer to kilometers. 
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The two observers will, however, disagree on the times of arrival of the 
satellite at various points in the coordinate system. Their clocks run at different 
rates, since they are in relative motion and in different parts of the gravitational 
field. But the time intervals recorded by the clocks will differ from each other 
by a factor of about GM £ /c 2 r, or about one part in 10 8 for the earth. Hence 
when the two observers compute the perigee shift, say, they will disagree on 
when the shift was a radians after a day's time by about a millisecond. This 
demonstrates that the timing problem is not important when tracking artificial 
satellites. 


ALTERNATIVE VIEWPOINT 

We can look at what we have done from the more usual non- Euclidian view- 
point. To do so requires a few words about coordinates. Let us confine our 
remarks to the (r, 0) plane, i.e. ® = 7T /2. 

Now r measures radial distance from the center of the earth; in fact, r = 
(proper area of sphere/ 4 tt) 1/2 (Misner, Thorne, and Wheeler, 1973; pg. 596). 

4> measures the angle on a sphere. So to get to the point (r,0) in physical space 
we go out the appropriate distance r and swing through the angle 0. 

We of course wish to find the set of values ( r, 0} which gives us the path of 
the satellite through physical space. We can accomplish this by doing the follow- 
ing: take a point (r, 0) in physical space and assign it a point (r E , 0 E ) in a 
Euclidian plane with numerical values r = r E and 0=0. Do this for all the 
points in physical space. Solve the equations of motion of the satellite by the 
methods of celestial mechanics to get its path in this Euclidian plane. In par- 
ticular, we wind up with a set of x and y values {x £ , y £ } for the track of the 
satellite. 

To find its p ath in phy sical space, what do we do? Take each (x £ , y £ ) in the 
set and put r, =/ + y 2 , 0 £ = Arctan y £ /x £ . Then set r = r £ and 0 = 0 £ . Go 

out distance r from the center of the earth and swing through angle 0 . Do tills 
for all the points in the set to get the track of the satellite through physical 
space. 

Now if we let the speed of light approach infinity, then the trajectory of 
the satellite becomes the unperturbed ellipse. We say that the difference between 
the perturbed and unperturbed positions is the displacement of the satellite. 
Computing the distances between the positions and getting so many centimeters 
displacement (for the case of the earth) may be done in the ordinary Euclidian 
sense, since the effect of the curvature of space is so small over such short 
distances that it may be neglected. 
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CONCLUSION 


The general relativistic correction to the position of an earth satellite with 
small orbital eccentricity has been shown to be about 17 cm per revolution. 

This effect is at present too small to be separated from other perturbing 
influences, such as radiation from the earth and atmospheric drag. However, 
improved knowledge of satellite perturbations and small atmospheric drag 
may allow the relativistic effect to be measured from observations of the 
proposed new geodynamic satellite LAGEOS. If it hoped that such observations 
will provide tests of the Einstein and Brans-Dicke theories of relativity. 
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APPENDIX 


This program takes the Keplerian elements of the perturbed and unperturbed 
orbits and calculates displacements. Explanations of its operation are given 
in the program itself. Some sample output is given. 


c 

C*****VAIN PROCfijlt' 

c 

C*****THIS PROGRAM CCMPUTES A SECUENCE OF DIFFERENCE VECTORS BETWEEN THE 
C POSITION IN AN UN PERT UP t)E 0 ORBIT AND THE POSITION IN THE ORBIT AS 

C FEFTOR06C EV THE RELATIVISTIC POTENTIAL AS T|-E MEAN ANCMALY OF 

C THE UNPEFTLREEC OnaiT INCREASES »ITH TIME. THE NLM0EF CF SUCH 

C VECTORS IS NP C l NT AND THE STEP S I ZE IN THE MEAN ANCMALY IS OEL . 

C THE COORDINATES OF THE DISTURBED POSITION RELATIVE TC THE 

C LND ISTUHFEC FCSITICN ARE LISTEC AND PLOTTED. THE SAME IS DONE FOR 

c THE VELOCITY VECTORS. EXCEPT THAT NO PLOT IS GIVEN. 

C THE KEPLERIAN ELEMENTS CF 0OTP PERTURBEC AND LNPERTUPEED ORBITS 

C ARE INITIALLY CHOSEN TO BE THE SAME AT PERIGEE. 

C 

C 

C*****NOTAT ICN FCR INPUT CATA. 


C 

INDEX TELLS WHAT OBJECT THE SATELLITE IS ORBITING AGCUT. INDEX=l 

C IS The SLN. 2 THE EARTH. AND 3 ANY U T HER GIVEN CEJECT. 

C ACM IS THE SEMIMAJOR AXIS CF THE ORBIT IN CENTIMETERS. 

El IS THE ECCENTRICITY. 

C XMASS3 IS THE M A 5 S OF THE OBJECT IN GRAMS IF INDEX=3. 

C D I S3 IS THE DISTANCE F AC T CR IN CENTIMETERS CHCSEN 5UCH THAT 

C ACM/D 15 3=1 APFPOX IMATELY , IF INDEX=3. 

NT E R M TELLS taHAT TERMS ARE TC EE RETAINEC IN THE PERTURBATION. 

C NTERM =0 GIVES THE FULL PER T UR E A T I CN . NTERM=l GIVES PERIODIC TERMS 

C CNLY « WHILE N T 6FM =2 GIVES SECULAR TERMS CNLY. 

C.....NPOINT IS THE NUMBER CF PCINTS IN THE SECUENCE. 

C XMSTRT IS THE INITIAL VALUE CF THE MEAN ANOMALY IN DECREES. 

C * . • . . DEL IS THE STEP SIZE OF THE MEAN ANOMALY IN DECREES. 

NPERI TELLS WHAT PGWEPS CF ECCENTRICITY ARE RETAINEC IN THE 

C ECCENTRICITY FUNCTIONS AND THEIR DERIVATIVES. IF NPEPT=1. THEN 

C POWERS OF ECCENTRICITY GREATER THAN 1 ARE NEGLECTEO. IF NPE«T=0» 

C THEN POWERS CHEATER THAN 2 ARE NEGLECTED. 

C 

C*****NOTAT I CN FCC CUTPUT CATA. 


C 

AND Y ARE THE CARTESIAN COPDINATES IN CENTIMETERS CF THE 

C PERTURBED POSITION RELATIVE TO THE UNP£RTU«e£0 PCSITICN (WHICH IS 

C AT X=0 AND Y=C.) 

C R = DSQR T ( X 2 ♦ Y**2) IS THE DISTANCE FROM THE ORIGIN TC (X.YI. 

C.....XDT AND YCT ARE THE CARTESIAN VELOCITY DIFFERENCES IN CM/5EC 
C BETWEEN THE PEPTUROEO AND UNFERTLRBEO VELCCITIFS. 

V = OSORT ( XCT * * 2 ♦ YDT**2) IS THE MAGNITUDE OF THE VELCCITY 

C DIFFERENCE. 

C 


000 I 
0002 
0 C 0 3 
0004 
0 C C 5 
0006 
0007 

0 C C 8 
CCOR 
OC 10 

0011 


IMPLICIT FEflL*e ( A— H . G—Z I 
C I MEN SIGN *1(4C0>.YJ (4001 
CATA RAD/C .C l 74532925 1994300/ 

CATA 0IGC/e.67C-8/.C/2.SSeCHO/ 

CATA XMASSl/l .SB7D33/.XMASS2/5.S76D27/.D I SI / l .456013/ 

CATA DIS2/6. 376 151 154D0/ 

CATA PI.P2.F3/SH SUN t 5H 6 AR T H. 5HO T HER/ 

C*+***REAO IN THE INPUT CATA. 

1 READ (5.2. END =12) I NDEX , A CM . E I . XM A S S 3 .D I S3 
READ (5.16) NTERM , NPOI N T , XMSTfi T .DEL . NPER T 

2 FORMAT I IS .CIS .5. FlO. 5. 2015. 5) 

16 FORMA T <2 I5.2F10.5. 15) 

C**+**cet THE VALUES OF THE I NCL I N A T ICN . LO NG I T UOE OF NCCE . AND ARGUMENT 
C CF PERIGEE TC ZERO. SO THAT THE ORBIT LIES IN THE X Y PLANE ANO 


JCH ir H ju 


CP AXl 


jf" f }i l 1J N r l h T u ix l C G£ □ IT LIE! 


WITH THE FCINT OF PERIGEE CN THE POSITIVE SICE OF THE AXIS. 


00 12 


FI 1=0. 0C0 

00 13 


C MEG 1 =C .CCC 

0014 


FNODEl = 0 .oco 


C ***** CHECK TC SEE 

0015 


IF (INDEX - ■ 

?C 16 

3 

XMASS= XM AS S 1 

00 17 


D 1 SF A C = C IS 1 

00 18 


P=P1 

0019 


GO TO 6 

0020 

4 

XMA SS = XN A S S2 

0 C 2 l 


C I SFiA C = C 1 SZ 

0 C 22 


F = P2 

CC23 


GO TO 6 

0 0 24 

5 

XMASS=X M ASS3 

0025 


C 1 SF.AC = C I S3 

0 C 26 


F=P3 

0027 

6 

CONTINUE 
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C C26 
OC 29 

0 C 30 

3 C 2 1 

C 0 3 2 
0C33 

0 3 34 

0C35 

0 0 3 6 
0037 
0(536 
OO 39 
0 04C 

004 1 
0C42 
OC 41 
0044 
0C45 
0046 
0 C 4 7 

0048 

0049 
C C 50 
CCS 1 
0 C52 

0 CS 3 
C C 5 4 
0055 

0 C56 
0 C 5 7 

o C 5 e 

0C59 
0 060 
0061 
0062 
0063 

CC64 
C 06 5 


0 C 6 6 
0 C 67 
0066 


0069 
OC 70 
0 07 1 
0 C 72 
OC7 3 


OC 74 


0 C 7 5 


OC 76 
OC 77 
C C 76 
0 C 79 

oc eo 

0 C 8 1 
0 082 
0 C ft 3 
0064 
0 C 6 5 

0 C 86 
CC67 

0068 


C *****CMC2 = 0 16 G * MASS/(5PEED OF LlGHT)**2 IN CCS UN 1 1 5 . 

GMC2 = BIC0’*XMASS/( C * * 2 ) 

SQ=(niGG*xMASS ) 70 I SF AC 

C *****CUNV IS 4 CONVERSION FAC T CH FOR THE VELOCITIES. 

CONV=DSCRT (SC ) 

4 1= ACM 70 ISF AC 

C * ** ** Hi? I TE OUT TEE EEAOINGS. 
write (6.7» 1 

7 FORMAT (1E1> 

WRITE (6,6) P .XMASS, DI SFAC 

6 FORMAT ( 7/7 . | 0 X . 1 6MD I STURE I NO BCJC Y = » 1 X , A 5 « I 0 X , SH M A S S » , 0 1 5 • 5 . ( X . 2HG 

1 M » 10X, 1 6FCISTANCE F ACTOR 1 . 1 X .0 15.5. 1 X.2HCIM 
WR I TE (6,9) ACE. El 

9- FORMAT ( / / / / / , ICX «2HA = » IX. DIE. 5. IX.2FCM, 2CX«2HE=, I X .FI0.5) 

WRITE (6.14) NFCINT.OEL 

14 FORMAT {/// ,iax,7FNPOINT=.lX,I5,20X,4FOEL=»IX,FlC.S.lX. 7HOEGREESI 

IF (NTERM - I ) 17, 1 p, 1 9 

17 WR 1 TE (6 » < C ) 

20 FORMAT ( 7 7 7 • ICX. 1 7HFULL PERTOH E A T I ON ) 

GC TO 23 

18 WR I TE (6.2 1) 

21 FORMAT ( /// . 1 CX. 1 9HPER (OC 1C TERES CNLV) 

GO TO 23 

19 WRITE (6 ,2 J ) 

22 FORMAT { 7 / 7 , ICX, 18HSECULAR TERMS CNLV) 

23 CONTINUE 

IF (NPEKT - 1) 24,26,26 


24 WRITE (6,25) 


2 5 

FORMAT' ( 7/7 , 1 C X . 7 7HTERMS LP TO AND INCLUDING 
1 IN THE ECCEMFIC1TY FUNCTIONS) 

SECCNO 

CROER 

RETAINED 


C-0 TO 2e 




26 

WRITE (6,27) 




27 

FORMAT ( 7/7 , l 0 X , 77HT ERMS LP TC ANO tNCLUCING 
1 IN THE ECCENTRICITY FUNCTIONS) 

FIRST 

ORDER 

RETAINED 

28 

CONTINUE 
WRITE (6,29) 




29 

FORMAT ( 1 C X , 2 1 F ANC THEIR DERIVATIVES) 
WRITE (6.7) 

WRI TE (6.50 




3 C 

FORMAT (77/, ICX) 
WRITE (6,13) 




1 3 

FORMAT ( 3X , 1EN , 4X , 1 2HM£AN ANCM AL Y , 7 X , 1 HX • I 4 X 

, 1 FT , 14 X 

. IHR , 

l 3X , 3HXOT 


1 * 1 2X. 3HYDT , 13X , 1HV ) 

WRITE (6,15) 

1 5 FORMAT ( 10 X .SH(OEGRfcES) . 7X , 4H( CM ) » I IX ,4H (CM > , 1 IX ,4H (CX ) , 9X,8H(CM/S 
IEC). 7X, 8F ( CM/S EC ) .7X.6H (CM/SEC),//) 

CW****DO THE ITER AT ICN. 

CO 10 N= 1 . NFC I NT 

XMEAN= ( N- 1 ) *0£L + XMSTRT 

XM l =X ME A N ♦ R A D 

,E1 .XMl.CWEGI ARE THE UNPEHT Lft EEC ELEMENTS. 

C4+***A2»E2,XV2 , C ME G £ ARE THE P ER T UR EEC ELEMENTS. 

C 4***4 C ALL DELTEL TC GET THE RELATIVISTIC CCRRECTICNS TO TEE ELEMENTS. 
CALL OELTEC ( N TER M , N PERT , GMC2 . AC M , E l , XM 1 , DEL A , CEL E , C EL M , OELG MG I 
A2r A 1 ♦ CELA/C1SFAC 
E2=E1 ♦ DELE 
XM2 = XM 1 + C EL M 

CMEG2=0MEG1 ♦ CELCMG 

C ***** CALL OETFV TO CONVERT THE KEPLERIAN ELEMENTS TO CARTESIAN 
C COORDINATES AND VELOCITIES. 

CALL 0ETRV(A1 ,E 1 <F1 I , XMl .CMEGl ,FNCOEl,XA,YA»ZA,XCTA.YCTA.ZOTA,«A,V 
1 A ) 

CALL OETRVIA2 ,E2,FI 1 »XM2 .CMEG2.FNCDE 1, XE .V6.2E ,XCT8« YC T9 ,20 T0 ,RB,V 
1 E > 

C****»C OMPUT E TFE DIFFERENCES IN COORDINATES AND VELOCITIES. 

X=( XB-XA)*C15FAC 

V= (YB— YA)*DISF AC 

2= ( 213- Z A )*C ISF AC 

R2 = I X* • 2 ) ♦ (Y**2) ♦ (Z**2) 

fi=CSGRT (R2 ) 

XDT= ( XDTE - X CT A 14CDNV 

YDT=(YDTE-YCT7)*CCNV 

ZDT=(ZDTB-ZCTA)*CONV 

V2= ( XDT ** 2 ) + I YD T *4 2 ) + (ZDT**21 

V = OStlHT ( V2 ) 

C *****PUT THE FEINTS IN THE ARRAY FOR FLCTI. 

XI (N) sX 
Y 1 (N) =Y 

C****+*RITE OLT TFE MEAN ANOMALY. DISTANCES. AND VELOCITIES. 

WRITE (6.11) N . XME AN.X.Y » fi « XD T »YCT .V 
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0 C 89 

10 

CONTINUE 


r*****CACL PLCT1 TO PLOT THE PC 

0 0 90 


CALC PLCTHXl *V1 •NPOINT) 

0091 


GO TO 1 

C 0 92 

12 

CONTINUE 

0093 

1 1 

FORMAT IlX.I4t2X.F10. 4. 3X 

0094 


STOP 

0095 


END 


I NTS. 
•6015.5) 
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COO 1 


1C 02 
000 3 
0004 

0 005 
0006 
3CC7 
0 00 8 
ccco 
0010 

001 1 


0012 
0 C 1 3 

0014 

0015 
0 0 16 
0 C 1 7 
00 18 

0019 

0020 
002 1 
0 022 
0023 
0 C24 
0C25 
0026 
0 C 2 7 
0 C 2 & 

0029 

0030 
0021 
0 032 
0 033 

0034 

0035 
0C36 
0 C3? 
C C3E 

0039 

0040 

0041 

0042 
004 3 
0 044 
0 045 
004 6 
0 C 4 7 


5 LflRQUT I NE CEL TEL I NT ER I* . NPE RT ♦ GNC2 ♦ A 0 . E C . X MQ , CE L A . D ELE . DELM . DEL CM 
1 0 1 
c 
c 

C44444THIS 5UBFOL T I NE FINDS TFE RELATIVISTIC CCRRECTICNS TC THE 
C KEPLERIAN ELEMENTS. WHICH ARE TO BE ACOEC TC TEE CNFERTUREED 

C ELEMENTS 1C GIVE THE ELEMENTS CF THE RELATIVISTIC ALL V PERTURBED 

C CRB IT* 

C 

C*****LN( J EW HJHtEC ELEMENTS 
C 

C44444A0 15 TEE SEE IRAJOR AXIS IN CENTIMETERS. 

C44444E0 IS TEE ECCENTRICITY. 

C44444XM0 15 TEE MEAN ANOMALY IN RACIANS. 

C 

C*****PERTUREEC ELEMENTS 

c 

C 4 4* 4 4 C ELA IS TEE CCRRECTICN TC THE SEN I N A JGR AXIS IN CENT I NE T ERS • 

C*»**#CELE IS Tee CCRRECTICN TC THE ECCENTRICITY. 

C444+4CELM IS TEE CCRRECTICN TO TEE MEAN ANCMALY IN RACIANS. 

C *****CELQMG IS TEE CORRECTION TO TEE ARGUMENT OF PERIGEE IN RADIANS. 
C44444THE PERTURE at ICNS IN I NCL I NAT l UN AND LONGITUDE GF NCCE ARE 
C ZERO ANC FENCE NOT INCLUDED IN TEE SUER CL T I NE . 

C44444NTERM TELLS MEAT TERMS TO RETAIN IN THE PERTURBATION. NTER M= 0 
C GIVES TEE FULL PE R TUR 0 A T I ON . MERM=| GIVES PERIODIC TERMS CNLY. 

C WHILE N TE R N = 2 GIVES SECULAR TERMS CNLY. 

C44444GMC2 = BIGG * MASS/ISPEEC OF LICFT 1**2 IN CGS UMTS. 

C44444NPERT TELLS WHAT POWERS OF ECCENTRICITY ARE RE T A l NE C IN THE 
C ECCENTRICITY FUNCTIONS AND THEIR DERIVATIVES. IF NPERT=1, THEN 

C POWERS OF ECCENTRICITY GREATER THAN I ARE NEGLECTED. IF NPER T = 0 • 

C THEN POWERS GREATER THAN 2 ARE NEGLECTEC. 

C 

IMPLIC I T RE AL + € ( A-H.O-Z > 

F 1=1 .ODO - ( EC**2 ) 

F 2 = DSQR T I F I ) 

SM = OSIN ( >MC ) 

CM=DCOS (XMC ) 

XM2 = 2 .ODC 4 XMC 
XM3=3.0DC AXMC 
SM2=DS l N ( >N2 ) 

CM2=DCOS ( X M 2 ) 

SM3=D51 N ( > N 3 ) 

C 4*»**G2 1 1 AND G2I2 ARE ECCENTRICITY FUNCTICNS. 

C*****GP211. GF 2 12. AND GP2I3 ARE CERIVAT I V E5 CF ECCENTRICITY FUNCTIONS. 
IF ( NPER T .EQ, 1) GO TO 4 
G2 1 1 = 1 .5CC* EC 

G21 2=19.000/4 .ODO ) 4 1 EO**2 > 

GP 21 1 = 1 .SCO ♦ (81 .000/1 e -CDO) *< t0**2 > 

GP2 1 2 = 4 .5CC4EC 

GP2 1 3= ( 159 .CCC/ 16 .0003 *( EC* *2 I 
CO TO 5 

4 G2 1 I = ( 1 .50 C MEC 
G2 12 = 0. CC C 

GP2 1 1= I .SOC 
GP2 12= ( 4 .5CC ) 4EO 
GP2 1 3 = 0 .OOC 

5 CONTINUE 

E 1 =G2 1 1 *CM 
e2 = GP2 I l ASM 
E 3 = G2 1 1 *SN 
C 1 =G2 I 2*C M 2 
C2=GP2 l 24SM2/2 .000 
C3=G2 I24SM2/Z .CDO 
F 2= GP2 1 2 * S M 3/3 .OOO 
FA = GMC2*(4 .00 0 )*F1 

C£LA=FA*BI - FA4G211 4 F A*C l - FA4G212 

FE = GMC2* < 2 , OOC )4(F1442)/(£04a0> 

CELE = FE *B 1 FE4G21I 4 FE4C1 - FE*G212 

FMS=GMC2* (2.0CC»/<F2*AO) 

X M S =F M S * X N 0 

F M I =GMC 2* (2.CCC)*(F1442)/(E0*AC ) 

F M2 =GMC 2* ( 12. 000) 4F 1/AO 

CELM=XMS - FM 14(02 4 C2 + P2 ) 4 FM24(B3 4 C2) 

FOMGS = GMC2*(3. ODO )/(F l* AO > 

F 0 MG 1 =G NC 2 4 (2 »CDO)AF1*F2/(EO*AO I 
DELOMG=FCNGS**MO 4 FOMGl*(e2 4 C2 + PZ ) 

IF INTERN - 1> 3,2.1 

1 DEL M= XM S 

OELOMG=FCNCS4XNO 
DEL A=0 . OOC 
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0 C 48 


06LE-0 .000 


C C4<5 


CO TO 3 

+ P2 ) ♦ 

0050 

2 

DELM— — F M 1 * ( 8 2 * C2 

0C51 


C£LGMG-FC*G 1* < E2 + 

C2 > P2 ) 

0 C 5 2 

3 

CONTINUE 


0G5 3 


return 


0 C 54 


END 



FM2*< 63 


C3> 
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0001 


SUBROUT IN£ CETRV(CEA. CEE. CE INC. CfcMA.CECMEC, CECAP, CEx.CEY.OEZ 
l CEXDT tOEVOl . CEZOT .OEHMAG.CEVMAG) 


0002 
0002 
COCA 
0005 
C CC6 
0007 

coce 

QCC9 
C C I 0 
OCl 1 
0012 
0013 
00 1 A 
CO l 5 
0016 


0010 
0 0 19 
0020 
0 02 1 
0022 
0C23 
0024 
0 C 25 
0C26 
G 027 
0028 

0029 

0030 
0 02 1 
0022 
0022 
0 0 34 
0035 
0 C 36 
0 C 37 
0 C30 
0039 
0 C 4 0 
0 C 4 1 
0042 
0 04 3 
0C44 
0045 
0C46 

0047 

0048 
C C 4 9 


c 

c 

C*****THIS SUeRCLTINE CONVERTS KEPL 6 R l AN ELEMENTS TC CARTESIAN 
C COORDINATES ANC VELOCITIES. 

C 

C ***** INPUT 

C 

C 4 * * ** C EA IS tee DIMENSIONLESS SEMIMAJUR AXIS [USUALLY CEA=1 
C APPROXIMATELY.) 

C 4 * 4 4 4 C LE IS TEE ECCENTRICITY. 

C****4CE I NC IS TEE INCLINATION IN RADIANS. 

C 4 4 4 4 4 0 EM A IS TEE MEAN ANOMALY IN RACIANS. 

C44444CECMEG IS THE ARGUMENT OF PERIGEE IN RADIANS. 

C444440ECAP IS TEE LCNGITUOE OF NUDE IN RADIANS. 

C 

C 44444C.UTPUT 
C 

c 44444CEX.UEY , AND CEZ ARE THE CARTESIAN X Y Z COORDINATES. TEE ORBITEO 
C OBJECT IS AT THE ORIGIN. 

C44444CEHMAG IS TEE CISTANCE FRCM THE ORIGIN TC THE SATELLITE. 

C 4 4 4* 4QEXDT « O E Y DT « AND CEZDT ARE THE VELOCITIES IN THE X.Y, AND Z 
C CIRECTICNS. 

C 4 44 44 CE VMA G IS TEE MAGNITUDE OF THE VELCCITY. 

C 44*4*1 HE DISTANCES AND THE VELOCITIES MOST EACH eE MULTIPLIED BY 
C CONVERSICN FACTORS IN THE MAIN PROGRAM TC CONVERT TO CGS UNITS. 

C IF THE CONVERSION FACTOR FCR TEE DISTANCES IS CISFAC (IN 

C CENTIMETERS J . THEN THE CONVERSION FACTOR FOR THE VELOCITIES IS 

C CONV = DSCRT{eiC G * MASS/ CISFAC) IN CGS UMTS. 

c 

IMPLICIT REAL *E < A-H, C-Z . * > 

T 40P 1=6 .20216* 2071 79586CC 
TAMEAN = CEMA 
IF (TAMEAN) 10.11*10 
10 TAMEAN = DMCD ( TAMEAN + T W CF I , T 4C P I ) 

ECCAI = TAMEAN + OEE 4CS I N ( T AMt AN ) + . 5DO 40E E 4* 2* C S I N { 2 . DO * T AM E AN ) 

CO 13 12 = 1 , 100 

CIFF = (CEE*CS INI ECCAI > - E CC A l + 1 A M t AN I / < l . C 0- C E E *CC C S ( ECC A l ) J 
ECCAZ = ECCAI ♦ O IFF 
SECCA2 = CSIN(ECCA2> 

CIFF = CA6S(ECCA2 - TAMEAN - CtE4SECCA2> 

IF (DIFF - 0 . 1C - 1 31 16,16.13 

13 ECCAI = ECCA2 

l*RITE (6 .2 JCEMA, CMA, DIFF 

2 FORMAT ( Sx .£ 3HNC CONVERGENCE IN KEPLERS ECUATICN - SUEROUTINE GETRV 


1 

/3D 1 6 .a ) 



STOP 


1 1 

6CCA2 = TAMEAN 



SECCA2= C 5 I N ( E CC A 2 ) 

I 6 

X = DCCS (ECC A 2 ) 

- OEE 


SP = 1.D0 - CEE**2 
£0 = D3GPT { SP ) 

Y = 5Q4SECCA2 
TRUEA = APCTAN(Y.X) 

C TR-UE A = CCCS(TRUEA) 

STRUEA * ESIN (TRUEA) 

CERMAG = C EA 4 Sf / ( 1 .DOPOEE *CTRUE A I 
X = OERMAG4CTFLEA 

V = OEPMAG4ST PLEA 
SN* = CSlMCfcOMEGI 
CS* = CCCS(CEOMEG) 

SNCPw = CS I N (CECAP) 

CSCPW = DCCS (CECAP) 

SNI - DS IN(CELNC) 

CSt. = CCCS(CEINC) 

AC - CStt4CSCF* - SNW*SNCPh*CSi 
EC = ,C SH4 £ NCR 4 + SNW*CSCP* *CS 1 
CC = SN44SM 

FC =— SN44CSCP1* - CSW*SNCPfc*C5I 

GC =-SM + SNCPI» +■ CSW4CSCPM4CS1 

K — C S 1*4 S N I 

CE X = AC*) + FC4Y 

CEY = £3C 4 X + CC*Y 

CEZ = C C 4 X ♦ K*Y 

xxD 1 = -STRUEA 

Y YD 1 - CEE + C TRUEA 

CEXDT- AC4XXC1 + FC* Y YD 1 

OE YDT = BC*XXC1 + GC4YY01 

CE 2D T = CC4XXC1 ♦ HC * Y YD 1 
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0 0 50 
005 1 

0052 

0053 

0054 

0055 
0 C 56 
0C57 


QEVMAG 
SGSMSQ 
FMULT 
CEXOT = 
GEYOT = 
CEZDT = 
RETURN 
END 


= DSCRT<2 .DO/OERMAG - I.DO/OEA) 

= CSCFT <OEXOT**2 ♦ CEYDT**2 ♦ CEZCT**2) 
= CE VM AG/SCSMSO 
CEXCT * F MULT 
CEYCT * F MULT 
OEZCT * F MULT 
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occi 


FUIMCTI GN A RCT AMS.C) 

0002 


IMPLICIT F£AL*6(A— H, 0— Z * * 1 

0003 


Y=S 

OOC4 


x=c 

0CQ5 


IF ,<X ) 1C € * 10C.I 16 

0006 

1 00 

IF { Y ) 1C2 , 104.1 06 

0007 

102 

ARCTAN=4 . 71 2366980 3 84689D0 

0CC8 


GO TO 1 38 

0 CC 9 

1 04 

ARCTAN=C ,0C0 

0010 


RETURN 

001 l 

106 

ARCTAN=l • £70796 3267948S6DC 

0012 


GO TO 138 

0013 

1 C 8 

IF < Y 1 110*112*114 

0014 

l 10 

ADD=3 . 1 4 159265 35 89793D0 

OG 15 


GO TO 124 

0016 

1 12 

ARCTAN=3. 14 159265358979300 

0017 


GO TO 138 

00 18 

1 14 

ADD=3. 14 1 S 9 26 5: 358979300 

0019 


GO TO 132 

0020 

116 

IF { Y ) lie. 120*122 

002 1 

1 18 

ADD-6. 233 1 85 3 C 7 l 7 9 5860 0 

0022 


GO TO 132 

0 C 2 3 

120 

ARCTAN^O .0C0 

0024 


GO TO 138 

0025 

122 

ADD^O.COO 

0026 

1 24 

IF (DAQS(Y>— CAES(X) ) 126.128.130 

0027 

126 

ARCT AN = C A TAN(Y/X) 4 AOD 

0028 


GO TO 138 

0029 

128 

ARCTAN= 0.7 85396 16 339 744 62 CO 4 ACC 

0 C 3 0 


GO TO 138 

0031 

130 

ARCTAN^l .570796 32679489600 - DAT AN ( X / Y ) 4 ADO 

0032 


GO TO 138 

0 C 3 3 

132 

IF <OABS(Y) -- CABStxn 126,134,136 

0034 

134 

ARCTAN=~C .765 3981 63397448 2DO 4 AOD 

0 C 35 


GO TO 138 

0036 

1 36 

ARCTAN--1 . 57C 79632679489600 - DATAN(X/Y> 4 ADC 

0 C 37 

138 

RETURN 

0038 


END 
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0001 


SUBROUT [N£ FLCT1I K.YiN) 


0002 

0003 

0004 
0 CO 5 


0CC6 

0CO7 

ocoe 
ocoo 
0010 
0011 
OC 1 2 
00 13 
0014 
OC 1 5 

oc ie 

00 17 

o o i a 

OC 19 
0020 
0021 
0022 
0023 
0 024 
0 025 
0026 
0 C27 

o c 2 e 

0029 
00 30 
0031 
0022 
0 033 
0 034 

0035 

0036 

0037 
0 0 38 
0 039 
0C40 
004 1 

0042 

0043 
0 C 4 4 
0045 
0 C 46 
004 7 
0 C 48 

0 04 9 
0 0 50 

0051 

0052 


0054 
C C55 
OC 56 
G C 57 
0058 
0 C 59 
0 C 60 

oce i 

0062 

0063 

0064 
0 C 65 
C C 66 
C C67 


C 

C*****1HIS SUBROUTINE PLOTS A GRAPH CONSISTING OF N POINTS IN ARRAYS 
C X AND Y . 

C The X AXIS IS HORIZONTAL AND THE Y AXIS IS VERTICAL. THE AXES ARE 

C PRINTEO AS DOTS. THE SCALE OF THE DIAGRAM IS ALWAYS ACJU5TED TO 

C FIT ONTC A SINGLE SHEET CF CCMFUTER PAPER. WITH 61 SPACES UP AND 

C 101 ACROSS. THE DISTORTION DUE TO ThE DIFFERENT SPACING IN THE X 

C AND Y DIRECTIONS IS REMCVEO, SO THAT A SCUARE LOCKS LIKE A SQUARE. 

C 

IMPLICIT REAL*e (A-H. C-Z) 

DIMENS ICN X (N 1 ,Y ( NJ 
D I MENS ICN AP < 1 C I . 61 > 

CATA BLANK .COT .D. Oh/ 1H * 1 F . • IF C . 1 HO/ 

C***** XSCAL6 TAKES CARE OF THE DIFFERENT SFACING IN THE VERTICAL AND 
C HORIZONTAL CIFCCTICNS. 

XSCALE=£.25CC/5.0CO 
YMl N=Y ( 1 ) 

YM AX^ Y ( I ) 

X M I N= X ( 1 ) 

XMAX=X ( I ) 

EC ALE 1 = 100. CDC 

ECALE 2 = 6 C .CCC 

NSC 2 = SC AL E2 + C . I DO 

CO 24 1=2. N 

C=Y(I) 

IF CU-YN AX ) 12.13,13 

13 YMAX=Q 
GO TO 14 

12 CONTINUE 

IF (O-YMIM IE. 14. 14 
15 YM1N=Q 

14 CONTINUE 
R = X( I ) 

IF (R-XMAX > 22 .23.23 

23 XMAX=R 
GC TO 24 

22 CONTINUE 

IF (R-XVIM 2E .24*24 

25 XM1N=R 

24 CONTINUE 

YLG TH= Y M A X- Y M IN 
XLGTH=XM AX-XMIN 
F AC = SC ALE l/( XLCTHWXSC AL E ) 

YF=YLGTh*FAC 

IF {YF - SCALES) 26,28.27 
27 FAC=SC ALE 2/YLGTH 

26 CONTINUE 
•RITE (6,41) 

WRITE I 6 » 7 C ) 

70 FORMAT ( / / / . 1 0 X . 3 3H I NFOP M AT l C N FROM SUBROUTINE PLGT1.///I 



WRI TE 

(6.71) 





7 1 

FORMAT 

( 10X, 3 7 HY MAX 

IS 

The maximum 

Y VALUE 

AT T A INED . > 


WRITE 

16.72) 





72 

FORMAT 

( I 0 X , 3 THY MI N 

I S 

THE V IN l MUM 

Y VALUE 

ATT A INEC . ) 


WH I TE 

(6,72) 





73 

FORMAT 

( 10 X . 1 9HYLGTH 

IS 

YMAX-YW IN. ) 




WR I TE 

(6.74) 





74 

FORMAT 

( ICX . 7CHXMAX . 

XMIN 

AND XLGTH 

ARE THE 

ANALCGCUS EXPRESSIONS 


ICR The 

X DIRECT ICN. ) 






WR I TE 

( 6 , t £ ) 






75 FORMAT (ICX.41HFAC ADJUSTS THE GRAPH TO THE PROPER SIZE.,///) 

WRITE (6.29) 

2 9 FORMAT I/// / ,5X.4HYMAX,12X,4HYMIN.I0X,4HXMAX. 12X.4FXMIN* 10 X, 5HYLGT 
1H.10X.5HX LG TH ,11X, 3HF AC , / ) 

WHITE (6,£C) V N AX , YM I N, XM A X , >M IN , YLGTH , XLC TH ,FAC 
60 FORMAT (e01S.£) 

CO £7 I- 1 *N 
57 X( 1 ) = X ( I ) * XSC ALE 
X M I N— XMINAXSCALE 
XMAX=XMAX*XSC ALE 
DO 29 I = 1 , 1C 1 
CO 29 J= 1 .6 1 

29 AP ( I , J ) —BLANK 

Jt = (-XRIM*FAC f I .500 
J2 = I - YM I N 1 * FAC F I .500 

IE (J1 .GT. C .AND. J 2 .GT. C) GC TG 20 
GO TO 31 

30 CONTINUE 

CO 32 I = 1 . 1 C 1 



30 



0C68 

32 

A P( I , J2 ) = OCT 

0 C 69 


CO 33 J = 1 * 6 1 

OC 70 

33 

AP< J1 , J )=CCT 

OC 71 

31 

CONTINUE 

0072 


CO 34 1 = 1 .Is 

OC 73 


Jl=(X(I) - XMIM*FAC + 1.500 

0 C 74 


J2 = CY(I) - VWirv)*FAC + 1.5CO 

0 C 75 


AP ( J 1 » J£ ) = C 

00 76 

34 

CONTINUE 

C 0 77 


WR l TE <6*41 ) 

OC 78 

4 1 

FORMAT ( 1 H 1 1 

0 C 79 


WRITE (6*76 ) 

0080 

76 

FORMAT <///,lCX) 

0081 


OO 35 J= 1 , € 1 

ocea 


J 3=6 1 - <J-1) 

0 C 83 


1»R I TE <6*401 ( AP( I * J 3 ) * 1=1*101) 

OC 84 

35 

CONTINUE 

0085 

40 

FORMAT < 1CX. 10 1A1 ) 

0 C 86 


RETURN 

o ce 7 


END 


SS**n 


as* 


31 



□ [STUBS ( NC EG DY= EARTH MASS” 0.5S7600 2S GM CISTANCE FACTOR^ 0.C37820 09 CM 

A= O.30AAOO 11 CM E= 0.055CO 

f>POINT = 73 0EL = 5-O0C00 DEGREES 

FELL PERTURiATION 

TERMS UF TO AND I NCLC C 1 NC SECCNO ORDER BETA! NED IN TEE ECCENTRICITY FUNCTIONS 
AND THEIR DERIVATIVES 


SAMPLE OUTPUT SAMPLE OUTPUT 
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MEAN ANOMALY X Y XOT YD Y V 

( OfcCREES) (CM < CM! t CM] C CM/ SEC » JCM/'SEC* 4 CM/ SEC) 



INFORMATION FROM SU E PC L T i N E 


PLOT 1 


YMflx IS THE MAXIMUM V VALLE ATTAINED. 

Y»IN IS THE MINIMUM V VALLE ATTAINED. 

YLGTH IS yMAx-VMIN. 

X MAX, AM IN ANC XLGTH APE TEL ANALOGOUS L XPH t S S I L N S FCF THE X DIRECTION. 

e ac adjusts thf graph tc tee proper size. 


YM AX TWIN 

0.167 7CD 02 -C.ee9'i7D 01 


X MAX 

C . ISZBoU 02 


XM l N 

-0. 1 74 2 2L 01 


YLGTH 

0.2SA6SC 02 


XLGTH 

0.170290 02 


FAC 

0.2J55eD Cl 
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